# Streamlined package management
pkgs <- c("dplyr", "broom", "haven", "tidyverse", "ggplot2")

# Function to load or install packages
usePackage <- function(p) {
  for (pkg in p) {
    if (!require(pkg, character.only = TRUE)) {
      install.packages(pkg, repos = "https://cloud.r-project.org/")
      library(pkg, character.only = TRUE)
    }
  }
}
# Load or install all required packages
usePackage(pkgs)
data <- read_dta("Propaganda_Experiment_working.dta") 
data$treatment <- ifelse(data$Treatment != "Control", 1, 0)
  


##############Main Outcome#############
# Define the variable names
H0_vars <- c("COVID_Policy_Support", "COVID_Preference_pc", "COVID_Protest", "Protest_Rightfulness")
control_vars <- c("")

# Conduct the linear regression analysis for the soft reopening and store the results in a list
coefs_list_soft <- list()
for (var in H0_vars) {
  reg_out_soft <- lm(as.formula(paste(var, "~ soft", paste(control_vars, collapse = " +"))), data = data, subset = (data$Treatment == "Soft" | data$Treatment == "Control"))
  reg_out_soft_soft <- lm(as.formula(paste(var, "~ soft_soft", paste(control_vars, collapse = " +"))), data = data, subset = (data$Treatment == "Control"| data$Treatment == "Soft_Soft"))
  coefs_list_soft[[var]] <- bind_rows(tidy(reg_out_soft), tidy(reg_out_soft_soft))
}

# Combine the coefficients for the soft reopening into a single data frame and filter for the "soft" and "soft_soft" terms
coefs_df_soft <- do.call(rbind, coefs_list_soft)
coefs_df_soft$variable <- factor(rep(names(coefs_list_soft), each = nrow(coefs_list_soft[[1]])))
coefs_df_filtered_soft <- coefs_df_soft[which(coefs_df_soft$term %in% c("soft", "soft_soft")), ]

# Conduct the linear regression analysis for the hard reopening and store the results in a list
coefs_list_hard <- list()
for (var in H0_vars) {
  reg_out_hard <- lm(as.formula(paste(var, "~ hard", paste(control_vars, collapse = " +"))), data = data, subset = (data$Treatment == "Hard" | data$Treatment == "Control"))
  reg_out_hard_hard <- lm(as.formula(paste(var, "~ hard_hard", paste(control_vars, collapse = " +"))), data = data, subset = (data$Treatment == "Control"| data$Treatment == "Hard_Hard"))
  coefs_list_hard[[var]] <- bind_rows(tidy(reg_out_hard), tidy(reg_out_hard_hard))
}

# Combine the coefficients for the hard reopening into a single data frame and filter for the "hard" and "hard_hard" terms
coefs_df_hard <- do.call(rbind, coefs_list_hard)
coefs_df_hard$variable <- factor(rep(names(coefs_list_hard), each = nrow(coefs_list_hard[[1]])))
coefs_df_filtered_hard <- coefs_df_hard[which(coefs_df_hard$term %in% c("hard", "hard_hard")), ]

# Combine the two data frames into a single data frame
coefs_df_combined <- rbind(coefs_df_filtered_soft, coefs_df_filtered_hard)
coefs_df_combined$propaganda_type <- factor(rep(c("Soft", "Hard"), each = nrow(coefs_df_filtered_soft)))
coefs_df_combined$group <- ifelse(coefs_df_combined$term %in% c("soft", "hard"), "Reopen", "ZCReopen")
coefs_df_combined$protest <- ifelse(grepl("Protest", coefs_df_combined$variable), 1, 0)
coefs_df_combined$varl = factor(coefs_df_combined$variable, labels = c("Assessment of Government Performance", "Prefer to reopen", "Willingness to Protest", "Protest Rightfulness"))

# FIGURE 2
ggplot(coefs_df_combined[coefs_df_combined$protest==0,], aes(x = variable, y = estimate, ymin = estimate - 1.96 * std.error, ymax = estimate + 1.96 * std.error, color = propaganda_type, linetype = group)) +
  geom_pointrange(alpha = 0.7, size = 0.2, position = position_dodge(width = 0.5)) +
  geom_hline(yintercept = 0) +
  coord_flip() +
  ylab("Coefficient Estimate of Propaganda") +
  theme_bw() +
  guides(color = guide_legend(title = "Propaganda Type"), linetype = guide_legend(title = "Policy Direction")) +
  scale_color_manual(values = c("red", "blue"), labels = c("Hard", "Soft")) +
  scale_linetype_manual(values = c("dotted", "solid"), labels = c("Reopen only", "Zero-COVID and Reopen")) +
  facet_wrap(~ varl, ncol = 1, scales = "free_y") +
  scale_x_discrete(labels = c(
    "COVID_Policy_Support" = "Soft Propaganda \n\n\nHard Propaganda", 
    "COVID_Preference_pc" = "Soft Propaganda \n\n\nHard Propaganda")) +
  theme(text = element_text(size = 15))
ggsave("Figure2.pdf", width = 10, height = 5)

# FIGURE 3
ggplot(coefs_df_combined[coefs_df_combined$protest==1,], aes(x = variable, y = estimate, ymin = estimate - 1.96 * std.error, ymax = estimate + 1.96 * std.error, color = propaganda_type, linetype = group)) +
  geom_pointrange(alpha = 0.7, size = 0.2, position = position_dodge(width = 0.5)) +
  geom_hline(yintercept = 0) +
  coord_flip() +
  ylab("Coefficient Estimate of Propaganda") +
  theme_bw() +
  guides(color = guide_legend(title = "Propaganda Type"), linetype = guide_legend(title = "Policy Direction")) +
  scale_color_manual(values = c("red", "blue"), labels = c("Hard", "Soft")) +
  scale_linetype_manual(values = c("dotted", "solid"), labels = c("Reopen only", "Zero-COVID and Reopen")) +
  facet_wrap(~ varl, ncol = 1, scales = "free_y", labeller = label_value) +
  scale_x_discrete(labels = c(
    "Protest_Rightfulness" = "Soft Propaganda \n\n\nHard Propaganda",
    "COVID_Protest" = "Soft Propaganda \n\n\nHard Propaganda")) +
  theme(text = element_text(size = 15))
ggsave("Figure3.pdf", width = 10, height = 5)

